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I.  Introduction 

When  solving  any  problem  in  fluid  mechanics  by  finite  difference  methods, 
the  first  order  of  business  is  to  generate  a  grid  which  fits  the  boundaries 
in  the  problem.  A  boundary  fitted  grid  is  essential  to  produce  an  accurate 
solution  because  the  region  in  the  vicinity  of  solid  surfaces  generally 
determines  the  character  of  the  flow  [11*. 

This  report  addresses  the  problem  of  generating  a  surface-fitted  grid  in 
a  pump  collector.  This  device,  sketched  in  Figure  1,  consists  of  an  annular 
inlet  ring  connected  to  a  doughnut-shaped  chamber  to  which  is  attached  a 
radial  exit  tube.  _ ^ 

\ 

Because  a  pump  collector  is  generally  a  three-dimens iohal  geometry,  the 
corresponding  surface-fitted  grid  is  also  three-dimensional.  ~hTo  develop  a 
capability  for  computing  such  a  grid  would  require  considerable  time  and 
effort.  The  development  time  can  be  reduced  drastically  by  making  two 
simplifying  assumptions  about  the  collector  geometry. v  The  resulting  geometry 
will  nevertheless  be  effective  for  conducting  parametric  studies  of  collector 
flows.  The  simplifying  assumptions  are  as  follows: 

(1)  The  inlet  ring  and  chamber  as  assumed  to  be  axially  symmetric. 

This  means  that  the  main  portion  of  the  device  is  formed  by 
rotating  the  "inlet-chamber  cross  section"  in  the  x-r  plane 
about  the  9  axis. 

(2)  In  the  r-9  plane,  the  cross  section  of  the  exit  tube  is 
assumed  to  be  rectangular  instead  of  round.  Also,  the 

*Numbers  in  brackets  designate  references  [see  p.  19] 
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profile  in  the  x-r  plane  is  assumed  constant.  Hence,  the 
pipe  is  formed  by  rotating  the  x-r  profile  through  an 
angle  A6. 

Thus,  a  two-dimensional  grid  can  be  fitted  to  the  cross  sections  of  the 
simplified  geometry  and  then  rotated  about  the  9  axis  to  produce  a  three- 

dimensional  grid. - * — '<*-  /  *-*  /j> 

Two  basic  approaches  are  available  for  generating  two-dimensional 
surface-fitted  grids.  These  are  differential  methods,  represented  by  the  work 
of  Thompson  and  co-workers  [2]  and  algebraic  methods,  represented  by  the  work 
of  Eiseman  [3] .  We  shall  extend  the  method  of  Amsden  and  Hirt  [4]  because  of 
its  simplicity  and  ease  of  application  to  the  present  problem.  In  the  Amsden- 
Hirt  procedure,  the  physical  coordinates  are  taken  to  be  solutions  in  the 
transformed  plane  of  a  linear  elliptic  system  which  consists  of  a  modification 
of  Laplace’s  equation  with  Dirichlet  boundary  conditions.  This  technique  is 
well  suited  to  our  ultimate  purpose  —  to  compute  frozen  vorticity  flows  in 
a  collector  —  because  it  generates  a  smoothly  varying  grid  in  the  physical 
coordinates  with  a  slowly  varying  mesh  cell  size. 

II.  Grid  Generation  Method 
Inlet-Chamber  Grid 

The  simplified  inlet-chamber  cross  section  is  formed  by  a  circle  and  a 
rectangle  with  one  side  of  the  rectangle  tangent  to  the  circle,  as  shown  in 
Figure  2.  A  simple  surface-fitted  grid  can  be  generated  for  this  cross 
section  by  transforming  the  boundary  into  two  connected  rectangles.  The 
original  cross  section  is  called  the  physical  plane  and  the  transformed 
rectangular  cross  section,  the  mapped  plane.  Inside  the  transformed  cross 
section  we  overlay  a  network  of  rectangular  grid  lines  of  uniform  spacing  h. 
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which  transforms  to  a  system  of  curved  grid  lines  in  the  physical  plane.  Thi 
process  of  generating  the  grid  can  be  broken  up  into  two  steps: 

(1)  We  prescribe  the  correspondence  between  points  on  the 

boundary  in  the  mapped  plane  with  those  in  the  physical 
plane.  This  prescription  is  arbitrary  and  can  be  used 

to  some  extent  to  control  distribution  of  grid  lines 

in  the  physical  plane. 

(2)  Next,  we  provide  a  prescription  for  determining  the 
location  of  grid  line  intersections  (vertices)  in  the 
physical  plane  which  correspond  to  vertices  in  the 
mapped  plane . 

Thus,  the  process  of  grid  generation  proceeds  backwards  from  the  mapped  plane 
with  a  rectangular  boundary  and  grid  to  the  physical  plane  with  a  curved 
boundary  and  grid.  We  will  now  cover  the  two  steps  of  the  process  in  detail. 

The  inlet-chamber  cross  section  and  grid  can  be  conveniently  described 
by  the  following  parameters: 

Tq  *  radius  of  chamber  circle  in  x-r  plane 

NLC  *  number  of  grid  intervals  on  chamber  radius 

Na  *  number  of  intervals  in  width  of  inlet 
*  number  of  intervals  in  length  of  inlet 
We  start  with  a  square  in  the  mapped  (£,n)  plane  with  origin  at  the 
square  center,  as  shown  in  Figure  3.  Since  the  dimensions  in  the  mapped 
plane  are  arbitrary,  for  convenience  we  choose  the  square  to  have  a  side 
of  length  two.  Then  the  step  size  h  in  the  mapped  plane  is 
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and  the  collector  inlet  length  l  and  vidth  w  in  the  mapped  plane  are 


A  ■  h  •  N_ 


w  ■  h  •  N 


A  ’ 


Using  an  (i,j)  index  system  to  designate  (5,n)  vertices,  we  have 


“  1  +  (i-1 )h  ,  1  <  i  <  NQ  +  1  , 


n  -  -  1  +  (j-l)h  ,  i  <  j  <  N  +  1  , 


where 


N  -  2Nt 


ng“n  +  nin  • 

On  the  boundary  in  the  mapped  plane,  corner  points  and  critical  control 
points  are  designated  by  the  letters  A-H  (see  Figure  3).  Critical  control 
points  in  this  case  are  E  and  F  which  .are,  respectively,  the  point  of  tangency 
of  the  circle  with  the  lower  inlet  segment  and  the  point  on  the  lower  inlet 
segment  directly  beneath  the  junction  of  the  circle  with  the  upper  inlet 
segment.  We  now  must  determine  corresponding  points  A-H  in  the  physical  plane 
as  well  as  all  intermediate  boundary  vertices  corresponding  to  intermediate 
boundary  vertices  in  the  mapped  plane. 

In  the  physical  plane,  we  choose  the  origin  of  coordinates  (x,r)  at  the 
center  of  the  circle.  The  location  of  vertices  on  the  bounding  circle  follows 
the  simple  prescription  of  Amsden  and  Hirt.  As  shown  in  Figure  4,  we  super¬ 
impose  the  square  with  side  of  length  two  on  a  circle  with  radius  r^  with 
common  origin  0.  Then,  a  straight  line  is  passed  through  the  origin  and  a 
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given  boundary  vertex  on  the  square  with  coordinates  (£,n).  .  The  intersection 

b 

of  this  line  with  the  circle  determines  the  corresponding  boundary  vertex 

(x,r)  in  the  physical  plane, 
b 

The  equation  of  the  line  is 


'  r 


where 


(6) 


a  *  (7)  =  given  , 

4  b 


and  the  equation  of  the  circle  is 


2  2  2 
x  +  r  =  rQ 


The  solution  of  these  two  equations  gives  (x,r)^: 


b  7T7T 

(1  +  a  ) 


(7) 


(8) 


(9) 


rb  ”  a  xb  ’  (10) 
where  the  sign  of  Xv  is  the  same  as  the  sign  of  .  If  £,  =  0,  then  we  have 

DO 

the  special  case 
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where  the  sign  of  r,  is  the  same  as  n. •  Thus,  all  of  the  vertices  in  the 

D  D 

physical  plane  between  A  and  E  are  determined  from  Eqs.  (9)  through  (12). 

Between  boundary  vertices  E-F-G-H-A,  associated  with  the  inlet  section, 
we  space  points  uniformly.  Thus,  between  E  and  F  the  boundary  vertex 
coordinates  are  given  by 

xb  “  5(f)  -  c  xA  ,  o  <  5  <  i  ,  (i: 

A 


Since  vertex  F  is  chosen  to  lie  directly  beneath  vertex  A  in  the  physical 
plane,  the  lengths  of  sides  FG  and  AH  are  the  same.  Hence,  for  uniform 
spacing  the  x  coordinates  of  vertices  on  these  boundary  segments  are  the 
same .  Thus , 

Xb  =  XF  +  ro("  ‘  *  XF  +  ro(C  “  l)  ’  1  <  ^  ^  ,  (1 


and 

r 

-  Tq  on  FG 

r .  c  n  AH 
v  A 

On  segment  GH,  the  coordinates  are 
xb  “  XG  ’ 


(1 


(1 


.ru  -  r. 


+  r„ 


*  ’  rH  +  <"  ■  •  "a  s  " > - '  • (1 
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With  the  inlet  corner  vertices  F-G-R-A  located  in  the  physical  plane, 
we  can  compute  the  inlet  dimensions.  The  inlet  length  s,  from  Figure  2, 
is 

s  =  xQ  -  xp  *  rQ(CG  -  l)  ,  (1 

and  from  Eqs.  (1)  and  (4) 


«G 


NIN  +  * 


'LC 


LC 


(21 


Hence,  the  dimensionless  inlet  length  in  terms  of  the  input  parameters  is 


s 


IN 


LC 


(2 


From  Figure  2,  the  inlet  width  t  is  seen  to  be 


r0  +  rA 


(2 


where  r^  is  obtained  from  Eq.  (10)  evaluated  at  boundary  vertex  A,  viz. 


rA  ’  aAXA 


aAr0 


(1  +  a‘) 
A 


(2 


From  Eqs.  (5)  and  (7),  we  find  that  is 


“A  *  nA  =  -  1  +  i 


(2 


LC 


Then,  combining  Eqs.  (23)  and  (24),  the  dimensionless  inlet  width  is 


-9- 


4  October  1984 
GHH: lhz 


t 


1 


Na  2 

[1  +  -  Y~) 

nlc 


1/2  * 

] 


Solving  Eq.  (25)  for  N^/N^g,  we  obtain 


1 


t/r  r 


[1  -  (1  -  t/r  )2] 


1/2 


(25) 


(26) 


The  inlet-chamber  geometry  dictates  that  NA  and  both  must  be  positive, 
i.e.,  Na/N^c  >  0.  When  this  ratio  is  zero,  we  find  that 

=  1  -  i-  =  0.2928932 
0  /2 

which  is  the  minimum  admissible  value  of  t/rg  for  this  mapping  scheme.  This 
lower  limit  is  a  consequence  of  how  the  vertices  on  the  square  are  mapped  to 
vertices  on  the  circle  and  corresponds  to  the  situation  when  line  OP  in 
Figure  4  intersects  the  lower  right-hand  corner  of  the  square.  This  limit 
can  be  removed  by  using  a  different  scheme  to  locate  the  vertices  on  the 
circle . 

In  the  determination  of  the  location  of  interior  vertices,  our  main 
interest  is  to  provide  a  prescription  that  will  give  a  smooth  variation  of 
vertex  locations  and,  consequently,  a  smooth  variation  of  the  metric 
coefficients.  For  the  present  application  to  the  calculation  of  inviscid 
flow  fields  where  flow  gradients  are  not  large,  the  simplest  such 
prescription  will  be  suitable. 
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Amsden  and  Hirt  chose  a  simple  prescription  for  locating  an  interior 
vertex  as  the  average  position  of  its  eight  nearest  neighbors.  A  still 
simpler  formula,  which  Amsden  and  Hirt  showed  gives  nearly  identical  results 
to  the  "eight  formula,”  is  to  place  the  interior  vertex  at  the  average 
location  of  its  north-south  and  east-west  neighbors.  Thus,  if  we  define  a 


two-component  column  vector  z  by 

*  9  J 


!i,j  =  [x’r]i,j 


(27) 


then  the  vertex  coordinates  by  the  latter  formula  are  given  by 


Ji,j  7  (*i-l,j  +  zi,j-l  +  zi+l,j  +  zi,j+l)  * 


(28) 


Equation  (28)  is  the  so-called  unit  square  formula  and  is  the  finite 
difference  solution  using  central  difference  approximations  of 


2  2 

8  z  ,  3  z  n 

— 2  +  ~  ’  0 

35  3n 


on  a  uniform  grid  in  5  and  n. 

In  the  location  of  interior  vertices,  we  will  use  the  unit  square 


formula.  Amsden  and  Hirt  solved  for  the  z  by  point  relaxation,  but  a  more 

*  J 


efficient  method  is  to  put  Eq.  (28)  into  tridiagonal  form  so  that  single  line 
overrelaxation  (SLOR)  may  be  used.  Thus,  along  i  =  constant  lines,  we  have 


Zi,j-1  "  4zi,j  +  Zi,j+1  3  Ri,j  * 


(29) 


where  R  is  a  two-component  column  vector  of  known  quantities  given  by 
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’  *1-1 ,J  "  Zi+l,j 


(30) 


In  the  relaxation  sweep,  for  the  direction  of  i  increasing,  the  latest  known 

values  of  z  are  used  to  evaluated  R  . 

1  f  J 

Along  line  i  specification  of  the  boundary  conditions  at  j  *  1  and 
j  *  jfflax  (boundary  vertices)  closes  the  tridiagonal  system.  A  solution  of 
this  linear  system  of  equations  is  then  obtained  by  the  usual  L-U 
decomposition  procedure  [5]. 

By  performing  numerical  experiments,  the  optimum  relaxation  factor  for 
the  collector  grid  was  found  to  be  1.5.  This  value  was  used  in  all  ensuing 
calculations.  Two  examples  are  presented  here  for  grids  with  step  sizes 
suitable  for  inviscid  calculations.  Table  I  lists  the  parameters  for  these 
examples . 


Case 

r0 

nlc 

nin 

na 

CO 

rt 

o 

t/r0 

1 

1.0 

10 

10 

4 

1.0 

0.4855 

2 

' 

1.0 

10 

10 

1 

10  ' 

1 

1.0 

1.0 

Table  I.  Collector  Grid  Geometry  Parameters. 


In  both  examples  the  number  of  iterations  required  for  convergence  was  28 
and  the  CPU  time  on  the  PSU  IBM  370/3033  was  2.3  seconds.  Convergence  was 
considered  achieved  when  the  maximum  change  in  x  and  r  during 

*  9  J  *  *  J 


succeeding  relaxation  sweeps  was  less  than  10 
Figure  5  and  Case  2  in  Figure  6. 


-4 


Case  1  is  shown  in 
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Exit  Pipe  Grid 

The  exit  pipe  profile  in  the  x-r  plane  is  assumed  to  be  symmetric  about 
x  *  0  and  is  made  up  of  three  straight-line  segments  plus  a  circular  arc,  as 
shown  in  Figure  7.  The  exit  pipe  boundary  A'-B'-C'-D'  is  transformed  into  a 
rectangle  over  which  is  laid  a  uniform  rectangular  mesh.  The  steps  for 
generating  this  grid  are  similar  to  those  for  generating  the  inlet-chamber 
grid. 

To  describe  the  exit  pipe  profile  and  grid,  the  following  parameters  are 

required  in  addition  to  r  and  N  : 

0  LC 

Nqw  *  number  of  intervals  in  half  width  of  exit  pipe  profile 

■  number  of  intervals  in  length  of  exit  pipe 

6  *  ratio  of  exit  pipe  final  width  #2  to  initial  width 

[see  Figure  7] . 

A  restriction  on  the  geometry  of  the  exit  pipe  profile  is  that  points  A' 
and  D'  must  lie  between  points  B  and  C  on  the  chamber  arc.  The  reason  for 
this  is  that  B  and  C  transform  into  the  corners  of  a  square  in  the  mapped 

plane  and,  hence,  A'  and  D'  cannot  lie  outside  these  points.  This  restriction 

means  that 

Now  <  \c  • 

The  half  width  of  the  exit  pipe  profile  at  the  chamber  arc  inter¬ 
section  is  given  by 


ui  ■  v 


(»  *  vz) 


2^71  • 


(31 ) 
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x  * 


(39) 


r  =*  r 


B' 


(40) 


The  interior  vertices  are  located  using  the  unit  square  formula,  Eq. 

(28),  and  the  resulting  tridiagonal  system  of  equations  is  solved  by  SLOR  with 
sweeps  in  the  £  direction. 

Two  examples  of  exit  pipe  profile  grids  are  presented  for  the  geometric 
and  grid  parameters  given  in  Table  II. 


Case 

r0 

nlc 

nol 

Now 

e 

1 

1.0 

- 

10 

10 

10 

1.0 

2 

1.0 

10 

10 

10 

0.5 

Table  II.  Exit  Pipe  Profile  Grid  Geometry  Parameters. 


In  both  cases,  the  number  of  iterations  required  for  convergence,  using  the 
same  criterion  as  for  the  collector,  was  17  and  the  CPU  time  was  0.9  seconds. 
Case  1  is  shown  in  Figure  8  and  Case  2  in  Figure  9. 

III.  Application  to  Calculation  of  a  Flow  Field 

In  this  section  we  shall  discuss  in  general  terms  how  the  previous  grid 
generation  scheme  can  be  used  in  the  numerical  calculation  of  a  collector  flow 
field.  We  start  with  the  strong  conservation  law  form  of  the  equations  of 
fluid  mechanics,  which  in  cylindrical  polar  coordinates  can  be  written  as 
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3E  3F  3G 
3x  3r  30 


H 


(41) 


where  E,  F,  G,  and  H  are  vectors  made  up  of  the  various  conserved  flow 
quantities.  New  independent  variables  5  and  n,  which  produce  a  grid  fitted  to 
the  inlet-chamber  cross  section  and  exit  pipe  profile  in  the  x-r  plane,  are 
introduced  by  the  transformation 

5  «  5(x,r)  ,  n  -  n(x,r)  .  (42) 


We  want  to  show  that  Eq.  (41)  can  be  transformed  into  a  strong  conservation 
form  in  the  (5»n,  9)  coordinate  system,  which  can  then  be  used  as  the  basis  for 
numerical  calculations. 

The  derivatives  with  respect  to  (x,r,9)  are  related  to  those  with  respect 
to  (5,n,9)  by  the  chain  rule  relations: 


3_ 

3x 


3 

nx  3n 


(43) 


3_  c  3  3 

3r  =  ^r  35  \  3n 


(44) 


where  subscripts  denote  partial  differentiation  with  respect  to  the  particular 
variable.  Upon  application  of  the  chain  rule  relations,  Eq.  (41)  becomes, 
after  rearrangement 


r.  3E  .  r  ae  \  .  r  at  .  (M  ,  * 
^x  35  ?r  35  ^nx  3n  nr  3n'  39 


3F 


3E 


3F ' 


3G 


H 


(45) 


For  later  convenience,  let  us  replace  5X»  etc*  to  x^*  xn>  etc.  by  use 
of  the  chain  rule.  If  we  operate  on  x  in  Eqs.  (43)  and  (44),  a  pair  of  linear 
equations  for  xF  and  x_  results  for  which  the  solution  is 
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Jx£  ’  nr  ’ 

(46) 

Jx«  “  “  5,  , 

(47) 

where  J  is  the  Jacobian  of  the  transformation  defined  by  the  following 
determinant: 


5x  \ 


£  n 

r  r 


5  n  ■  C  n 
x  r  r  x 


(48) 


Similarly,  by  operating  on  r  in  Eqs.  (43)  and  (44),  we  get  a  pair  of  linear 

equations  for  r  and  r  which  has  the  following  solution: 

£  n 


Jcc  ■  -  \  • 

(49) 

Jr  -  5 

(50) 

Equations  (46)  through  (50)  are  used  to  replace  £x*  etc.  by  x^,  etc.  in 
Eq.  (45).  Thus,  Eq.  (45)  becomes 

3F 


jfr  —  -  x  —  1  +  jf-  r  —  +  x  —1  +  —  -  H 

1  n  3£  n  3£;  1  £  3n  £  3n;  30 


To  simplify  Eq.  (51),  the  following  identities  are  needed: 


3E  3  r  r.  'i  _ 
r  T7  *  -r;  !r  E  -  Er. 


-  x 


3F 
n  3£ 


3£ 


>nF)  +  'V, 


-  r 


3E 
£  3n 


3n  *■  £ 


(r.E)  +  Er 


£n 


3F  3  f  _  \  _ 

x£  3n  "  3n  ^X£F'  "  Fx 


£n 


,***»•»  ^  •  *  «  •  »  *  «  '  *  *  ■  •  »•  «*•*•*•  ,  •  •  «  .  •*■%** 


.*.**.«.  ■  .  ^  ■  A. 


(51) 


(52) 
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Then,  Eq.  (51)  becomes 

3  3  3G 

J  85  ('nE  -  V)  +  J  (-  rsE  +  +  ae  •  H  •  <”> 

The  final  transformed  strong  conservation  form  is  obtained  from  Eq.  (53)  by 
dividing  through  by  J,  noting  that  J  ■  J(C.n)  and  defining  the  following 
transformed  vectors: 


E  - 


r_E  -  x_F 

n  n 


(54) 


F  *  -  r^E  +  x^F  , 


(55) 


G  -  GJ 


(56) 


H  -  HJ 

Then,  the  result  is 


»  +  if  +  J£  „  H 

35  3n  30  H 


(57) 


(58) 


Equation  (58)  thus  holds  in  a  coordinate  system  for  which  the  solid  boundaries 

are  families  of  the  coordinates.  Also,  from  Eq.  (58)  an  appropriate  finite 

difference  algorithm  can  be  derived  which  will  be  described  in  a  later  report. 

In  the  numerical  solution  of  Eq.  (58),  the  metric  coefficients  x_,  x  , 

5  n 

r^  and  r^  are  required  in  Eqs.  (54)  through  (57).  Since  the  step  size  in  the 
C-n  plane  is  a  constant,  simple  central  difference  expressions  of  second-order 
accuracy  can  be  used  to  compute  these  quantities,  viz. 
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OJ  “  +  0(h2)  ,  (59) 

*  i.j  Z 

(*)  -  +  0(h2)  ,  (60) 

n  ltj  2h 

where  z  is  defined  by  Eq.  (27).  These  derivatives  need  to  be  computed  only 
once  for  a  given  geometry  and  then  stored  for  later  use  in  the  fluid 
mechanical  calculation. 
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Fig.  1.  Geometry  for  a  Pump  Collector. 


Collector  Grid,  Case  1 
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Fig.  7.  Exit  Pipe  Profile  Geometry  and  Grid. 
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